
library("rstudioapi")     
##Set working directory
setwd(dirname(getActiveDocumentContext()$path))
#setwd("../")
#setwd("../")
getwd() 
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(xlsx)
library(estimatr)

# Load final Data file: 
load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

# added by AGB: centre and scale age and income
library(dplyr)
final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

#########################################
## lm.cluster {miceadds}
#########################################
# Model L_7b with Ncash_cat variable
L_7b <- glm.cluster(vaccine_intention ~ Ncash_cat +
                      dist2 + dist3 + dist4 + dist5 + dist6, # fixed effect for district
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) #### individual treatment

# Model L_7c with cash_cat, Gender, Education level and Age
L_7c <- glm.cluster(vaccine_intention ~ cash_cat + Gender + Agec + 
                      Education + Incomec + Employed + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
L_7c2 <- glm.cluster(vaccine_intention ~ individual_treatment + Gender + Agec + 
                       Education + Incomec + Employed + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalV2, family=binomial,
                     cluster = final_finalV2$Q123) 

summary(L_7c)
exp(cbind("Odds ratio" = coef(L_7c), confint.default(L_7c, level = 0.95)))

# Model L_7c with cash_cat, Gender, Education level and Age ## its the same as 7C! 
# gets dropped and goes to the other table 
L_7d <- glm.cluster(vaccine_intention ~ Ncash_cat + Gender + Agec + 
                      Education + Incomec + Employed + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
summary(L_7d)
exp(cbind("Odds ratio" = coef(L_7d), confint.default(L_7d, level = 0.95)))

# Model L_8c with cash_cat, Gender, Education level and Age
L_8c <- glm.cluster(vaccine_reported_combo ~ cash_cat  + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec +
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
L_8c2 <- glm.cluster(vaccine_reported_combo ~ individual_treatment  + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec +
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalV2, family=binomial,
                     cluster = final_finalV2$Q123) 
summary(L_8c)
exp(cbind("Odds ratio" = coef(L_8c), confint.default(L_8c, level = 0.95)))

# Model L_8d with cash_cat, Gender, Education level and Age
L_8d <- glm.cluster(vaccine_reported_combo ~ Ncash_cat + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
L_8d2 <- glm.cluster(vaccine_reported_combo ~ individual_treatment + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalV2, family=binomial,
                     cluster = final_finalV2$Q123) 
summary(L_8d)
exp(cbind("Odds ratio" = coef(L_8d), confint.default(L_8d, level = 0.95)))

# Models 9 with the new data: 
L_9c <- glm.cluster(ActVacApril ~ cash_cat  + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123)
L_9c2 <- glm.cluster(ActVacApril ~ individual_treatment  + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalV2, family=binomial,
                     cluster = final_finalV2$Q123)
summary(L_9c)
exp(cbind("Odds ratio" = coef(L_9c), confint.default(L_9c, level = 0.95)))

L_9d <- glm.cluster(ActVacApril  ~ Ncash_cat + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) #### individual treatment
summary(L_9d)
exp(cbind("Odds ratio" = coef(L_9d), confint.default(L_9d, level = 0.95)))

## Export results for model L_7d and L_7c into a Stargazer table: 
tab_7c1 = exp(cbind("Odds ratio" = coef(L_7c), confint.default(L_7c, level = 0.95)))
tab_7c2 = exp(cbind("Odds ratio" = coef(L_7c2), confint.default(L_7c2, level = 0.95)))
tab_7d1 = exp(cbind("Odds ratio" = coef(L_7d), confint.default(L_7d, level = 0.95)))
tab_8c1 = exp(cbind("Odds ratio" = coef(L_8c), confint.default(L_8c, level = 0.95)))
tab_8c2 = exp(cbind("Odds ratio" = coef(L_8c2), confint.default(L_8c2, level = 0.95)))
tab_8d1 = exp(cbind("Odds ratio" = coef(L_8d), confint.default(L_8d, level = 0.95)))
tab_9c1 = exp(cbind("Odds ratio" = coef(L_9c), confint.default(L_9c, level = 0.95)))
tab_9c2 = exp(cbind("Odds ratio" = coef(L_9c2), confint.default(L_9c2, level = 0.95)))
tab_9d1 = exp(cbind("Odds ratio" = coef(L_9d), confint.default(L_9d, level = 0.95)))

rownames(tab_7c1) = c("Intercept", "Cash", "Health", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_7c2) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_7d1) = c("Intercept", "Cash", "Cash_Placebo", "Health", "Health_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8c1) = c("Intercept", "Cash", "Health", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8c2) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8d1) = c("Intercept", "Cash", "Cash_Placebo", "Health", "Health_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_9c1) = c("Intercept", "Cash", "Health", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_9c2) = c("Intercept", "Health", "High Cash", "Low Cash", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_9d1) = c("Intercept", "Cash", "Cash_Placebo", "Health", "Health_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")

# Convert function, to create outputs for stargazer tables: 
convert_fun = function(data, name_mod){
  new_vec = c()
  for (i in c(1:dim(data)[1])){
    new_vec[i] = as.character(paste(format(round(data[i,1],2), nsmall = 2),paste("(",paste(format(round(data[i,2],2), nsmall = 2), format(round(data[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
  }
  new_vec = data.frame(new_vec)
  rownames(new_vec) = rownames(data)
  colnames(new_vec) = name_mod
  return(new_vec)
}

col_7c1 = convert_fun(tab_7c1, "Model1")
col_7c2 = convert_fun(tab_7c2, "Model2")
col_7d1 = convert_fun(tab_7d1, "Model3")
col_8c1 = convert_fun(tab_8c1, "Model1")
col_8c2 = convert_fun(tab_8c2, "Model2")
col_8d1 = convert_fun(tab_8d1, "Model3")
col_9c1 = convert_fun(tab_9c1, "Model1")
col_9c2 = convert_fun(tab_9c2, "Model2")
col_9d1 = convert_fun(tab_9d1, "Model3")

tb1 = merge(col_7c1, col_7c2, by = 0, all = TRUE);rownames(tb1) = tb1[,1];tb1[,1] = NULL
tb2 = merge(col_8c1, col_8c2, by = 0, all = TRUE);rownames(tb2) = tb2[,1];tb2[,1] = NULL
tb3 = merge(col_9c1, col_9c2, by = 0, all = TRUE);rownames(tb3) = tb3[,1];tb3[,1] = NULL
ttb1 = merge(tb1, tb2, by = 0, all = TRUE);rownames(ttb1) = ttb1[,1];ttb1[,1] = NULL
ttb2 = merge(ttb1, tb3, by = 0, all = TRUE);rownames(ttb2) = ttb2[,1];ttb2[,1] = NULL
colnames(ttb2) = c("Intention", "Intention V2", "Reported", "Reported V2", "Actual", "Actual V2")
no.obs = c((L_7c$glm_res$df.null),(L_7c2$glm_res$df.null), (L_8c$glm_res$df.null), (L_8c2$glm_res$df.null), (L_9c$glm_res$df.null), (L_9c2$glm_res$df.null))
aic_mod = round(c(L_7c$glm_res$aic, L_7c2$glm_res$aic, L_8c$glm_res$aic, L_8c2$glm_res$aic, L_9c$glm_res$aic, L_9c2$glm_res$aic),0)
log_lik = round(c(as.vector(logLik(L_7c$glm_res)), as.vector(logLik(L_7c2$glm_res)), as.vector(logLik(L_8c$glm_res)), as.vector(logLik(L_8c2$glm_res)), as.vector(logLik(L_9c$glm_res)), as.vector(logLik(L_9c2$glm_res))),0)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(ttb2)

# Manuscript: Table 3: Odds Ratio from Regression of Vaccine Outcomes on Treatments and Covariates 
tab2_manuscript = rbind(ttb2, mod_spec)
tab2_manuscript = tab2_manuscript[c("Cash", "High Cash", "Low Cash","Health", "Access (+5)" ,"Age (+10 yrs)", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food (+50)","Social Media (+10)","District 2", "District 3", "District 4", "District 5","District 6", "Intercept", "Observations", "Akaike Inf. Crit.", "Log Likelihood"),]
stargazer(tab2_manuscript, type = "text", summary = FALSE)
stargazer(tab2_manuscript, type = "latex", summary = FALSE)

# shorter table, with simple "Yes" for Covaraites and Fixed effects
short_tab = tab2_manuscript
short_tab[5,] = rep("Yes", 6)
rownames(short_tab)[5] = "Covariates"
short_tab[14,] = rep("Yes", 6)
rownames(short_tab)[14] = "Fixed Effects"
short_tab = short_tab[-c(6:13, 15:19),]
stargazer(short_tab, type = "text", summary = FALSE)
stargazer(short_tab, type = "latex", summary = FALSE)


require(ggplot2)
tab_7c1 = as.data.frame(tab_7c1);tab_7c1$Variable = rownames(tab_7c1)
tab_7c1$Variable = factor(tab_7c1$Variable, levels = rev(c("Intercept", "Cash", "Health","Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed")))
tab_7c2 = as.data.frame(tab_7c2);tab_7c2$Variable = rownames(tab_7c2)
tab_7c2$Variable = factor(tab_7c2$Variable, levels = rev(c("Intercept", "Health", "High Cash", "Low Cash","Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed")))
tab_8c1 = as.data.frame(tab_8c1);tab_8c1$Variable = rownames(tab_8c1)
tab_8c1$Variable = factor(tab_8c1$Variable, levels = rev(c("Intercept", "Cash", "Health","Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)")))
tab_8c2 = as.data.frame(tab_8c2);tab_8c2$Variable = rownames(tab_8c2)
tab_8c2$Variable = factor(tab_8c2$Variable, levels = rev(c("Intercept", "Health", "High Cash", "Low Cash","Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)")))
tab_9c1 = as.data.frame(tab_9c1); tab_9c1$Variable = rownames(tab_9c1)
tab_9c1$Variable = factor(tab_9c1$Variable, levels = rev(c("Intercept", "Cash", "Health", "Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)")))
tab_9c2 = as.data.frame(tab_9c2); tab_9c2$Variable = rownames(tab_9c2)
tab_9c2$Variable = factor(tab_9c2$Variable, levels = rev(c("Intercept", "Health", "High Cash", "Low Cash", "Access (+5)", "Male", "Age (+10 yrs)", "High-Educate", "Medium-Educate", "Low-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)")))

p7c1 = ggplot(tab_7c1[-c(1,12:16),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 1: Vaccine Intention", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p7c2 = ggplot(tab_7c2[-c(1,13:17),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 2: Vaccine Intention Low-High", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p7c1
p7c2

p8c1 = ggplot(tab_8c1[-c(1,13:17),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 3: Reported Vaccine", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p8c2 = ggplot(tab_8c2[-c(1,14:18),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 4: Reported Vaccine Low-High", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p8c1
p8c2

p9c1 = ggplot(tab_9c1[-c(1,13:17),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 5: Actual Vaccine ", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p9c2 = ggplot(tab_9c2[-c(1,14:18),], aes(Variable, `Odds ratio`))+
  geom_point(size=3)+ theme(legend.position="none") + labs(title="Model 6: Actual Vaccine Low-High", x="", y="Odds Ratio and 95% CI")+
  geom_errorbar(aes(ymin=`2.5 %`,ymax=`97.5 %`),width=0.1) +  geom_hline(yintercept = 1, color="red") + coord_flip()
p9c1
p9c2

# Manuscript: Figure 2: Financial Incentive and Vaccination Outcomes: Odds Ratio
require(cowplot)
plot_grid(p7c1, p7c2,
          p8c1, p8c2,
          p9c1, p9c2,
          labels = NULL,  label_size = 12,  ncol = 2,  nrow = 3)


pdf("Main_Figure2.pdf", width = 10.0, height =9.5)
plot_grid(p7c1,p7c2,
          p8c1, p8c2,
          p9c1, p9c2,
          labels = NULL,  label_size = 12,  ncol = 2,  nrow = 3)
dev.off()

#### SPILLOVER DATA: 

final_finalSP = mutate(final_finalSP, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

# The Spillover models: 
L_7cs <- glm.cluster(vaccine_reported ~ Vcash_cat + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 
summary(L_7cs)
exp(cbind("Odds ratio" = coef(L_7cs), confint.default(L_7cs, level = 0.95)))

# Model 8cs: 
L_8cs <- glm.cluster(ActVacApril ~ Vcash_cat + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 
# Further tables for models:
tab_7cs = exp(cbind("Odds ratio" = coef(L_7cs), confint.default(L_7cs, level = 0.95)))
tab_8cs = exp(cbind("Odds ratio" = coef(L_8cs), confint.default(L_8cs, level = 0.95)))
rownames(tab_7cs) = c("Intercept", "VCash_cash", "VCash_Health", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)", "District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8cs) = c("Intercept", "VCash_cash", "VCash_Health", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)", "District 2", "District 3", "District 4", "District 5","District 6")

# Create proper output for stargazer: 
col_7cs = convert_fun(tab_7cs, "ModelSP1")
col_8cs = convert_fun(tab_8cs, "ModelSP2")

tsp1 = merge(col_7d1, col_8d1, by = 0, all = TRUE);rownames(tsp1) = tsp1[,1];tsp1[,1] = NULL
tsp2 = merge(col_7cs, col_8cs, by = 0, all = TRUE);rownames(tsp2) = tsp2[,1];tsp2[,1] = NULL
tbsp1 = merge(tsp1, col_9d1, by = 0, all = TRUE) ;rownames(tbsp1) = tbsp1[,1];tbsp1[,1] = NULL
tbsp2 = merge(tbsp1, tsp2, by = 0, all = TRUE) ;rownames(tbsp2) = tbsp2[,1];tbsp2[,1] = NULL
colnames(tbsp2) = c("Vaccine Intention", "Vaccine Reported", "Actual Vaccine", "Vaccine Reported SP", "Actual Vaccine SP")

no.obs = c((L_7d$glm_res$df.null), (L_8d$glm_res$df.null), (L_9d$glm_res$df.null), (L_7cs$glm_res$df.null), (L_8cs$glm_res$df.null))
aic_mod = round(c(L_7d$glm_res$aic, L_8d$glm_res$aic, L_9d$glm_res$aic, L_7cs$glm_res$aic, L_8cs$glm_res$aic),0)
log_lik = round(c(as.vector(logLik(L_7d$glm_res)), as.vector(logLik(L_8d$glm_res)), as.vector(logLik(L_9d$glm_res)), as.vector(logLik(L_7cs$glm_res)), as.vector(logLik(L_8cs$glm_res))),0)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(tbsp2)


tabSP_manuscript = rbind(tbsp2, mod_spec)
tabSP_manuscript = tabSP_manuscript[c("Cash","Health",  "Cash_Placebo" , "Health_Placebo", "VCash_cash", "VCash_Health" , 
                                      "Access (+5)", "Age (+10 yrs)", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food (+50)",
                                      "Social Media (+10)","District 2", "District 3", "District 4", "District 5","District 6", "Intercept", "Observations", "Akaike Inf. Crit.", "Log Likelihood"),]

stargazer(tabSP_manuscript, type = "text", summary = FALSE)
stargazer(tabSP_manuscript, type = "latex", summary = FALSE)

stargazer(tabSP_manuscript[,-1], type = "text", summary = FALSE)
stargazer(tabSP_manuscript[,-1], type = "latex", summary = FALSE)


short_tab = tabSP_manuscript
short_tab[7,] = rep("Yes", 5)
rownames(short_tab)[7] = "Covariates"
short_tab[16,] = rep("Yes",5)
rownames(short_tab)[16] = "Fixed Effects"
short_tab = short_tab[-c(8:15, 17:21),]
stargazer(short_tab, type = "text", summary = FALSE)
stargazer(short_tab, type = "latex", summary = FALSE)

###### Version with Low and High Cash instead of simple Cash and Cash Placebo: 
L_7d <- glm.cluster(vaccine_intention ~ Ncash_catV2 + Gender + Agec + 
                      Education + Incomec + Employed + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
tab_7d1 = exp(cbind("Odds ratio" = coef(L_7d), confint.default(L_7d, level = 0.95)))


L_8d <- glm.cluster(vaccine_reported_combo ~ Ncash_catV2 + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
tab_8d1 = exp(cbind("Odds ratio" = coef(L_8d), confint.default(L_8d, level = 0.95)))


L_9d <- glm.cluster(ActVacApril  ~ Ncash_catV2 + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) #### individual treatment
tab_9d1 = exp(cbind("Odds ratio" = coef(L_9d), confint.default(L_9d, level = 0.95)))

rownames(tab_7d1) = c("Intercept", "Health", "Health_Placebo", "HighCash", "HighCash_Placebo", "LowCash", "LowCash_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8d1) = c("Intercept", "Health", "Health_Placebo", "HighCash", "HighCash_Placebo", "LowCash", "LowCash_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_9d1) = c("Intercept", "Health", "Health_Placebo", "HighCash", "HighCash_Placebo", "LowCash", "LowCash_Placebo", "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)","District 2", "District 3", "District 4", "District 5","District 6")

col_7d1 = convert_fun(tab_7d1, "Model3")
col_8d1 = convert_fun(tab_8d1, "Model3")
col_9d1 = convert_fun(tab_9d1, "Model3")

# The Spillover models: 
L_7cs <- glm.cluster(vaccine_reported ~ Vcash_catV2 + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 
summary(L_7cs)
exp(cbind("Odds ratio" = coef(L_7cs), confint.default(L_7cs, level = 0.95)))

# Model 8cs: 
L_8cs <- glm.cluster(ActVacApril ~ Vcash_catV2 + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 
# Further tables for models:
tab_7cs = exp(cbind("Odds ratio" = coef(L_7cs), confint.default(L_7cs, level = 0.95)))
tab_8cs = exp(cbind("Odds ratio" = coef(L_8cs), confint.default(L_8cs, level = 0.95)))

rownames(tab_7cs) = c("Intercept", "VCash_Health","VCash_HighCash", "VCash_LowCash" , "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)", "District 2", "District 3", "District 4", "District 5","District 6")
rownames(tab_8cs) = c("Intercept", "VCash_Health","VCash_HighCash", "VCash_LowCash" , "Male", "Age (+10 yrs)", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food (+50)", "Employed", "Social Media (+10)", "Access (+5)", "District 2", "District 3", "District 4", "District 5","District 6")

# Create proper output for stargazer: 
col_7cs = convert_fun(tab_7cs, "ModelSP1")
col_8cs = convert_fun(tab_8cs, "ModelSP2")

tsp1 = merge(col_7d1, col_8d1, by = 0, all = TRUE);rownames(tsp1) = tsp1[,1];tsp1[,1] = NULL
tsp2 = merge(col_7cs, col_8cs, by = 0, all = TRUE);rownames(tsp2) = tsp2[,1];tsp2[,1] = NULL
tbsp1 = merge(tsp1, col_9d1, by = 0, all = TRUE) ;rownames(tbsp1) = tbsp1[,1];tbsp1[,1] = NULL
tbsp2 = merge(tbsp1, tsp2, by = 0, all = TRUE) ;rownames(tbsp2) = tbsp2[,1];tbsp2[,1] = NULL
colnames(tbsp2) = c("Vaccine Intention", "Vaccine Reported", "Actual Vaccine", "Vaccine Reported SP", "Actual Vaccine SP")

no.obs = c((L_7d$glm_res$df.null), (L_8d$glm_res$df.null), (L_9d$glm_res$df.null), (L_7cs$glm_res$df.null), (L_8cs$glm_res$df.null))
aic_mod = round(c(L_7d$glm_res$aic, L_8d$glm_res$aic, L_9d$glm_res$aic, L_7cs$glm_res$aic, L_8cs$glm_res$aic),0)
log_lik = round(c(as.vector(logLik(L_7d$glm_res)), as.vector(logLik(L_8d$glm_res)), as.vector(logLik(L_9d$glm_res)), as.vector(logLik(L_7cs$glm_res)), as.vector(logLik(L_8cs$glm_res))),0)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(tbsp2)


tabSP_manuscript = rbind(tbsp2, mod_spec)
tabSP_manuscript = tabSP_manuscript[c("Health", "HighCash","LowCash", "Health_Placebo", "HighCash_Placebo", "LowCash_Placebo", "VCash_Health", "VCash_HighCash", "VCash_LowCash" , 
                                      "Access (+5)", "Age (+10 yrs)", "Male", "High-Educate", "Medium-Educate", "Low-Educate","Employed","Mean-Food (+50)",
                                      "Social Media (+10)","District 2", "District 3", "District 4", "District 5","District 6", "Intercept", "Observations", "Akaike Inf. Crit.", "Log Likelihood"),]

stargazer(tabSP_manuscript, type = "text", summary = FALSE)
stargazer(tabSP_manuscript, type = "latex", summary = FALSE)

# Table 04: 
stargazer(tabSP_manuscript[,-1], type = "text", summary = FALSE)
stargazer(tabSP_manuscript[,-1], type = "latex", summary = FALSE)

short_tab = tabSP_manuscript
short_tab[10,] = rep("Yes", 5)
rownames(short_tab)[10] = "Covariates"
short_tab[19,] = rep("Yes",5)
rownames(short_tab)[19] = "Fixed Effects"
short_tab = short_tab[-c(11:18, 20:23),]
stargazer(short_tab, type = "text", summary = FALSE)
stargazer(short_tab, type = "latex", summary = FALSE)


stargazer(short_tab[,-1], type = "text", summary = FALSE)
stargazer(short_tab[,-1], type = "latex", summary = FALSE)


